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I. INTRODUCTION 

Quantum correlations in many-body and optical sys- 
tems are playing a crucial role in a variety of fields, from 
the microscopic study of novel states of matter in quan- 
tum fluids^, to the control and suppression of noise in 
optical systems 3 - - — , to the exploration of analog models 
of gravitational and cosmological systems^—. 

In this perspective, non-linear optical systems exhibit 
a rich variety of phenomen a 11 ' 12 , involving the interplay 
of parametric scattering with losses, non-linear energy 
shifts, and the peculiar dispersion of light in confined ge- 
ometries. As a most significant example of this physics, 
an example of phase transition of the Bose-Einstein class 
has recently started being investigated in the optical 
parametric oscillation in planar geometrie o 13 ' 14 . 

Theoretical work has suggested that further insight in 
this physics can be obtained from the spatial and tem- 
poral pattern of correlations of the emitted ligh t 15 ' 16 : as 
the critical point for parametric oscillation is approached, 
the correlation length and time of the parametric emis- 
sion show a divergence that is closely related to the one 
of thermodynamical phase transitions. 

To this purpose, semiconductor microcavities in the 
strong coupling regime appear as most favorable sys- 
tem o 13 ' 14 ' 17 as they are intrinsically grown with a pla- 
nar geometry, and nonlinear interactions between the 
dressed photons - the so-called polaritons - are remark- 
ably strong. First experiments demonstrating quantum 
correlations have recently appeared 1 ^—. 

In the present work, we investigate the physics of the 
intensity correlations that are generated by parametric 
scattering processes. Depending on the pump frequency 
and intensity, different regimes can be identified. In ad- 
dition to the usual short-distance correlations that are 
generally present in any interacting system, spontaneous 
parametric emission processes are responsible for addi- 
tional long-distance ones, whose correlation length di- 
verges as the optical parametric oscillation threshold is 
approached. In addition to their intrinsic interest, an ex- 
perimental study of intensity correlation in this simplest 
system will open the way to the investigation of more 



complex geometries that in many aspects mimic the be- 
havior of quantum fields in curved space-times^ 1 .. 

Our approach is based on a generalized, non- 
equilibrium Bogolubov-de Gennes approximation in 
which weak fluctuations around the classical coher- 
ent field are described in terms of a quadratic 
Hamiltonia n 14 ! 22 " — . Solving in frequency space the cor- 
responding quantum Langevin equations allows to obtain 
predictions for the spatial and temporal dependence of 
the in-cavity intensity correlations. These then directly 
transfer to the near-field correlations of the emitted light. 
The study of the simplest planar geometry is extended 
to the more realistic case of weakly disordered systems: 
in this regime, the effect of disorder is shown not to qual- 
itatively modify the peculiar correlation pattern. 

In Section HU we review the theoretical formalism 
based on the Bogolubov approximation and we discuss 
how correlations transfer from the in-cavity field to the 
emitted radiation. In Section IIIII A we present the nu- 
merical results for a one-dimensional, spatially uniform 
case. Approximate analytical calculations are presented 
in Section [HU B and used to physically interpret the nu- 
merical results. Numerical results for a two-dimensional 
disordered system are discussed in Section [TV] Conclu- 
sions are drawn in Section fVl 



II. FORMALISM 



In this section we briefly review the formalism which 
will be used to calculate the intensity correlations of the 
emitted light. As we are restricting to the case of small 
fluctuations around a strong coherent field, a Bogolubov 
approach is adapted to describe the dynamics of quantum 
fluctuations. In the simplest case of a spatially homoge- 
neous geometry the Bogolubov equations can be worked 
out to analytical expressions for the physical observables. 
In the general case, numerical results can be obtained by 
inverting the Bogolubov matrix for the specific geometry 
under investigation. 
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A. The system Hamiltonian 

We describe the exciton and the photon quantum fields 
^(r) and ^(r) as scalar Bose fields. The Hamiltonian 
is the sum of terms describing the free propagation of 
excitons and photons Hq, their dipole coupling Hr, the 
two-body exciton-exciton interaction Hx, saturation of 
the exciton-photon coupling H s and the external pump- 
ing Hj&2k 

H = H +H R + H x + H s +H p . (1) 
The free propagation term has the form 

H = f d 2 v¥ x {v)[E x - {h 2 /2m x )V 2 +U x (r)]$ x (r) 

d 2 r^(r)[e c (-iV) + C/ c (r)]^ c (r) (2) 



where E x is the exciton energy, m x its effective mass, and 
e c is the photon dispersion in the planar cavity; U x , c are 
the external potentials acting on respectively the photon 
and the exciton. The dipole exciton-photon coupling has 
the form 



H R = m R / d 2 r[*t( r )* c (r)+h.c] 



(3) 



and is quantified by the Rabi frequency Qr. The effective 
two-body exciton-exciton interaction term 

H x = \ 9x J rf 2 r^(r)*t (r)^(r)#,(r) (4) 

models both Coulomb interaction and the effect of Pauli 
exclusion on electrons and hole a 26 ' 27 . 

H s = g s J d 2 r[i>i(r)& x (r)y x (r)y x (r) + h.c] (5) 

is the term modeling the saturation of the exciton oscil- 
lator strength 2 ^. The external pump term has the form 



H p = hj d 2 r[*t(r) F p {r, t) + h.c. 



(6) 



In the following of the paper, we assume that the sys- 
tem is driven by a continuous- wave monochromatic pump 
with a plane wave spatial profile 



F(r,t) 



-iwptgikp-r ^ 



(7) 



Excitons and photons decay in time with a rate ^ x and 
7c respectively. 



B. Bogolubov-de Gennes formalism 

To make the problem analytically tractable, we per- 
form the Bogolubov approximation 2 ^: the two quantum 



fields are split into a strong coherent -classical- compo- 
nent and weak quantum fluctuations 



*x(c)(M) = e 



$o;(c)(r) + # x ( c )(r,t) 



(8) 



This decomposition is then inserted into the Hamiltonian 
and all terms of third and higher order in the fluctu- 
ations are neglected. This leads to a quadratic Hamilto- 
nian for the fluctuation field that can be attacked with 
available theoretical tools. 

The classical component $> x ( c )(r) is obtained from the 
Heisenberg equations of motion of the quantum field 
i& x / c )(r) by factorizing out the multi-operator averages. 
This leads to the following pair of generalized Gross- 
Pitaevskii equations: 



fkJ p <b x (r) 



ft 2 V 2 

■- + U x {v) - ih lx /2 + g x | $ x (r) | 2 

2m x 

2 ffs Re{$;(r)$ c (r)}) 

(M} R + g s I $ K (r) | 2 ) $ c (r), (9) 
e c (-iV) + U c (r) - ih lc /2~\ * c (r) (10) 
Ml R + g s I $ a (r) | 2 ) * x (r) + hF e ik - T , 



Note that the classical field § C}X (r) that is obtained 
from this equation does not include the correction due 
to the backaction of fluctuations onto the coherent com- 
ponent^. To take this effect into account, one should in- 
clude the contribution (5ip x (r) Sip x (r)) of the fluctuating 
fields to the density and then iteratively solve (|9lflO|) up 
to convergence. As in the present paper we are interested 
in the correlation properties of the field fluctuations, this 
effect can be safely neglected in what follows. 

Within the input-output formalism^ 2 -, the quantum dy- 
namics of fluctuations can be written in terms of quan- 
tum Langevin equations of the for m 22 ' 23 ' 29 



where 



ihd t S^{r,t) = MS^(r,t) + M(r,t), (11) 



8* = WxMiMcMtF (12) 



is the four-component quantum fluctuation field. The 
matrix M is obtained by linearizing the Heisenberg equa- 
tion of motion for the quantum field around the classical 
component. In our case, it has the form: 



T — i ^ 7x 
J-x I 2 



M 



MR 

V -VIZ 



X X \ * 

12 ) 



\^xx 
^12 




x ' 2 


-(£-)* 


^12 




-Ob 








-T c — i 



h . 



(13) 
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with 



Tx(r) 



h 2 V 2 

Irrir, 



+ U x (r) -hw p + 2g x | $ x (r) 





+ 4 9s Re{$*(r)$ c (r)}, 


(14) 


t(r) 


= e c (-iV) + C/ C (r) , 


(15) 


Eff(r) 


= fl x *»(r) + 2fl.* !B (r)* c (r) ) 


(16) 


n«(r) 


= Ml R + 2g s | ^(r) | 2 , 


(17) 




= S^(r). 


(18) 



The quantum Langevin force 

f = (/x 1 /x! fci f c ) 



describes the zero-point fluctuations in the input field. 
As the field dynamics takes place in a small frequency 
window around uj p and the input field is assumed to be 
in the vacuum state, the spectrum of quantum Langevin 
force can be approximated by a white noise in both space 
and time, 



(f 6 (r,t)rt(r',t'))=~, x S 6 , x S{r-r')d(t-t'). 



(20) 



This approximation is accurate in the present case where 
no other radiation is incident onto the cavity in addition 
to the driving at u> p and the light-matter coupling MIr is 
small as compared to the natural oscillation frequencies 
E x and eJ2»2i. 



and making use of ([2~Tj) , simple expressions for the two- 
operator averages are obtained 

4 r 

g x (r,r',u J ) = J2 / M*> WV. 1>*<", W^r, 1)^ , 
i=i •* 

(27) 

4 r 

g c (r,r',Lo) = J2 / M*> 3>*<s, l\M^ l \v, 3>r, , 

i=i 

(28) 

(19) m x (r,r»=^ / ds{s, l\M~y, 2)*(s, l\M-y, l)r, , 

(29) 

4 /• 

^ / ds(s,?|Mj 1 |r',4)*(s,/|AC 1 |r,3)r i . 

(30) 



m c (r,r',w) 



i=i 



where T = 2irh 2 (0,-f x ,0,-f c ) T and the matrix 

M u = M -hull (31) 

is evaluated in the basis \(r,j)). Here, |r) spans the real 
space and the label j = 1, ...,4 refers to the block form 
of the Bogolubov matrix (fl"3")h 



C. Evaluation of the correlation functions 

In the present case of a monochromatic pump, tempo- 
ral homogeneity guarantees that the correlation functions 
only depend on the time difference t — t' and different fre- 
quency components are decoupled: 

hoj 5*(r, w) = M<5*(r, u) + M{r, u) . (21) 

The two-operator averages of the fluctuation fields can 
be written as 

~g^ c) (r,t;r',t') = J due^-^g^y^) , (22) 
and 

m« c) (M;r',t') = j 'dwe-^-^mg^r.r'.w). (23) 
where we have defined 
~g$> c) (r,r',u>) = J dJ {8^ x(c) {v , u)6^ x{c) (t', U ')) , (24) 

™g ) (ry iW )=/m ( .)(^)%)(«'X)>. (25) 

By introducing the frequency-domain correlation func- 
tions 

(/ 5 (r, w)/t(r', w '))f = 27r 7x ^, x 5(r - r^w - «') (26) 



D. Intensity correlations 

The present paper is focussed on the spatio-temporal 
correlations of the intensity fluctuations of the in-cavity 
photon field. As usual in generic planar microcavities, 
the in-cavity intensity correlations directly reflect into 
the near-field intensity correlation of the emitted light 
outside the cavity. An explicit calculation of the relation 
between the in-cavity and the emitted fields in the spe- 
cific case of semiconductor DBR microcavities is given in 
the Appendix. 



The four operator correlation function describing the 
correlation of the intensity fluctuations of the in-cavity 
photon field is defined as usual as: 

G^(r,t;r',t') = t)^, 0*c(r', t / )* c (r, t)) . 

(32) 

At the level of the Bogolubov approximation, three- and 
four-operator averages of the fluctuation field can be ne- 
glected and the correlation function (|32fl can be written 
in terms of two-operator averages as 

G^(r,r',t) = |$ c (r)| 2 |$ c (r')| 2 + |$ c (r')| 2 (^(r)4(r)> 

+ \Mr)\ 2 (^ c (r')Mr')) 

+ 2Rc{$ c (r')*$ c (r)(^(r ; i)4( r ' ;0 ))} 

+ 2Rc{$ c (r')*$ c (r)*(4(r,t)V; c (r',0))}33) 




FIG. 1: Real (a-g) and imaginary parts (b-h) of the Bogolubov dispersion as a function of the wavevector k, for the four regimes 
A-D marked in panel (i). The exciton interaction energy is g \$ x \ 2 = 1(T 4 meV (a, b), g\$ x \ 2 = 0.04 meV (c, d), g\$ x \ 2 = 0.08 
meV (e, f) and gl^l 2 = 0.61 meV (g, h). (i) Bistable loop in the in-cavity exciton field intensity as a function of the incident 
pump intensity. The red dashed line represents the unstable solutions. The points labeled with the letters A-D mark the four 
regimes considered in panels (a-h). 



A particularly useful quantity is the reduced correlation 
function: 

G^Hr,t-,r\t')= G(2 )^ T '^ -l, (34) 
n c (r) n c (r'J 

scaled by the intensity of light 

n c (r) = |$ c (r)| 2 + (^t (r) ^ c(r))5 (35) 

which can be written in the following simple form: 

G( 2 '(r,r',t) ~ 2Rc{$ c (r')*$ c (r)(V;t( r ,t)4( r ',0)) 

5 (r / )*$ c (r)*(^ c (r > t)^ c (r / > 0)>} X 



4V 



x [|$ c (r)| 2 |$ c (r')| 2 ]- ] 



(36) 



In the next section, we will study the spatial and tempo- 
ral features of this quantity in the most relevant cases. 



III. ID UNIFORM SYSTEM 
A. Numerical results 

In this section we investigate the case of a spa- 
tially uniform system under a plane wave monochromatic 
pump. Spatial homogeneity guarantees that the correla- 
tion functions only depend on the relative spatial coor- 
dinate r = x — x' . 

We use typical parameters for a GaAs microcavity 
with N = 10 quantum wells. In particular, we choose 
hflji = 10 meV, hr) c = fvy x = fry = 0.1 meV, and we take 
zero exciton-photon detuning at k = 0, E x = e c (k = 0). 
For the sake of simplicity, in the present section we fo- 
cus our attention on a quasi-one-dimensional geometry 
where polaritons arc transversally confined to a length 



Iid ~ 1 M m - This corresponds to effective ID nonlin- 
ear coupling constants g x = 1.5 x 10 _3 meV /im and 
g s = 0.5 x 10~ 3 meV ^m. Such geometries are presently 
under active experimental investigation^ 2 -. Extension to 
the 2D case does not qualitatively affect the physics and 
will be discussed in Sec. IIVI in connection to disorder 
issues. 

We consider a configuration where the pump wave vec- 
tor k p = 2.2 //m _1 is close to the magic wave vector, and 
the pump energy is Hu) p — e c (0) = —8.7 meV. This config- 
uration allows to study the most significant regimes by 
simply varying the pump intensity. In Fig. Q] (i) , we plot 
the corresponding exciton field intensity as a function of 
the pump intensity |-Fb| 2 : the hysteresis cycle typical of 
bistable systems is apparen t 25 i 33 i 34 . 

In the following, we will consider four different density 
regimes, marked on the bistability loop by the points 
A-D. Two of them (A, B) correspond to low-density 
conditions far from any instability, the point C is very 
close to the OPO threshold, and point D is in the stable 
region above the bistability and parametric oscillation 
region. 

In the present uniform system, the elementary excita- 
tions on top of the steady-state of the pumped system 
can be classified in terms of their wavevector k. Exam- 
ples of their Bogolubov dispersion 2 ^ are shown in Fig. Q] 
(a-h) for the four different regimes marked in Fig. [T] (i) 
by the points A-D. In the figure, we restrict the field of 
view to the lower-polariton region that is involved in the 
physics under investigation here. 

For increasing intensities, we observe: i) the low- 
density parametric regime A (panels a, b), where the Bo- 
golubov modes reduce to the single-particle dispersion; ii) 
the regime B corresponding to moderate densities (pan- 
els c, d), where the imaginary parts are modified and a 
small region of flattened dispersion appears at the cross- 
ing of the Bogolubov modes; iii) the regime C close to 
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OPO threshold (panels c, f), where the imaginary part 
of one mode tends to zero and a flat region in the Bo- 
golubov dispersion is apparent; iv) the non-parametric 
configuration D, for intensities larger than the bistabil- 
ity threshold, where the normal and the ghost dispersions 
are well separated and the eigenmodes of the system tend 
again to the single-particle dispersion, yet blue-shifted by 
interactions (panels g, h). 

For each of this A-D regimes, we have computed the 
spatio-temporal pattern of intensity correlations with the 
model described in the previous section. The results are 
presented in the next subsections. An analytical inter- 
pretation will be given later on in Sec. IIIII B. 



1. Low-density parametric luminescence 

In Fig. [5] (a), we show the spatio-temporal pattern 
of the intensity correlation in the low-density regime A 
corresponding to Fig. Q] (a, b): this is characterized by 
a system of parallel fringes and a butterfly-shaped enve- 
lope. The fringe amplitude is almost vanishing inside the 
cone delimited by the group velocities of the signal and 
the idler modes (represented by the thin solid lines in 
the figure) and largest on the edges of the cone. Further 
away in space and time, it decays back to zero as a con- 
sequence of the finite polariton lifetime. At zero delay 
t = 0, the central bunching peak around r = shows a 
weak narrow dip [Fig[2Jb)]. At larger delays [FigJSJc)], 
the correlation signal is weaker in between the dot-dashed 
lines indicating the edges of the cone. 



2. Moderate- density parametric luminescence 

In Fig. [3] (a), we show the spatio-temporal pattern 
in the moderate density regime B corresponding to the 
energy dispersion shown in Fig. Q] (c, d). As compared 
to the low-density case, the qualitative shape of the cor- 
relation pattern is qualitatively modified: the system of 
parallel fringes extends in a significant way into the in- 
terior of the butterfly shape and the exponential decay 
in the external region takes place at a slower rate. This 
latter effect is a direct consequence of the increased life- 
time of the Bogolubov modes in the parametric region 
(see Fig.ffl(d)). 



3. OPO critical region 

In Fig. @] (a) , we show the spatio-temporal pattern for 
the regime C very close to the OPO threshold. In this 
case, the system of parallel fringes extends to the whole 
(r, t) space and correlations are non- vanishing even at 
very long time and space separations. This is a conse- 
quence of the diverging correlation length of fluctuations 
in the critical region^. The zero-delay t = cut of the 
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FIG. 2: (a) Pattern of the reduced intensity correlations 
G^ 2 ' (r, t) as a function of the spatial distance r and of the 
temporal delay t, for the low-density regime A, corresponding 
to Fig. Q](a, b). The two thin lines, r = v s g t and r = vlt, high- 
light the group velocities of the signal and idler modes, (b) 
Cut of the figure (a) along the zero delay t = line, (c) Cut 
along the t — 5 ps line. In panels (b) and (c), the analytical 
prediction Eq. (|45[) is shown as a dashed line. In panel (c), 
the two vertical dash-dotted lines delimit the region where 
the analytical result (|45[) predicts vanishing correlations. 
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FIG. 3: (a) Pattern of the intensity correlations G™(r, t), 
for the moderate density regime B, still well below the OPO 
threshold and corresponding to Fig. Q] (c, d). The two thin 
lines, r = Vgt, and r = v g t, highlight the group velocities of 
the signal and idler modes. 
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FIG. 4: (a) Pattern of the intensity correlations & 2 \r,t) 
for the density regime C just below the OPO threshold and 
corresponding to Fig. [JJ(e, f). (b) Cut of the figure (a) along 
the zero delay t = line, (c) Cut along the t = 10 ps line. 
In panels (b) and (c), the dashed line is the result of the 
analytical prediction Eqs. (|47|) and (|48)l . 



correlation pattern is shown in Fig. |4] (b) and is charac- 
terized by a system of fringes centered at r = with a 
central bunching peak. The weak dip that was visible in 
the weak intensity regime is no longer present. 



4- Non-parametric regime 

For even larger values of the pump intensity beyond 
the bistability threshold, the system is pushed into an 
higher density D regime, where the blue-shift induced 
by interactions brings parametric scattering processes far 
off-resonance [Fig. [1] (g, h)]. As a consequence, the sys- 
tem of parallel fringes due to parametric correlations al- 
most completely disappears, as it can be seen in Fig. [5] 

(a) . The remaining correlation pattern has an opposite 
shape, with substantial correlations limited to the region 
inside the cone delimited by the maximal group velocity 
of the polariton dispersion. At zero delay t = [Fig. [5] 

(b) ] no spatial fringe is visible and the only feature is an 
anti-bunching dip at r = 0. At finite time delays [Fig. [5] 

(c) ], some spatial fringes appear, but they do not show 
any dominating periodicity. This suggests that a con- 
tinuum of processes at different fc's are simultaneously 
taking place, each of them characterized by a different 
phase velocity 
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FIG. 5: (a) Pattern of the intensity correlations G*- 2 ' (r, t) 
for the high-density regime D, corresponding to Fig. [1] (g, 
h). In this regime, the large blue-shift of the polariton 
modes prevents the parametric oscillation. The two thin lines 
r = ±vT ax t mark the maximal group velocities, (b) Cut of 
the previous figures along the zero delay t = line, (c) Cut 
along the t — 5 ps line. In panels (b) and (c), the analyti- 
cal prediction Eq. (|52|l (red dashed line) and Eq. (|53jl (green 
dash-dotted line) are also shown. In panel (c), the two vertical 
thin dash-dotted lines delimit the region where both the nu- 
merical result and the analytical result (152[) predict vanishing 
correlations. 

B. Analytic model and interpretation 

To physically understand the numerical results pre- 
sented in the previous section, approximate analytic for- 
mulas for the behavior of G' 2 ' can be extracted for the 
most significant limiting cases. To simplify the analysis, 
we restrict our attention to the lower polariton branch. 
This approximation is accurate as long as the Rabi split- 
ting is much larger than all other energy scales, e.g. 
Ml R > 7, MIr > g x \<Z> x \ 2 , g s \® c $x\, and is safely ful- 
filled in all the cases that we are considering here. In this 
regime, the dynamics of the system can be described by a 
single polariton field. Its dispersion, linewidth 7 and non- 
linear interaction coefficient g are immediately obtained 
from the linear eigenmodes of the Gross-Pitacvskii equa- 
tions (HJ and (fT0|) once interactions are neglected^. At 
the same level of approximation, also the fc-dependence 
of the photon Hopfield factor Ck — C can be neglected. 
As a result, the same Hopfield coefficient |C| 4 appears in 
both the numerator and in the denominator of the inten- 
sity correlation function (|34p of the photon field, which 
then reduces to the corresponding quantity for the lower- 
polariton field. 

For a spatially homogeneous system under a 
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monochromatic plane wave pump, equations can be writ- 
ten in the momentum space and the reduced second-order 
correlation function reads 



& 2 \v,t)=Gf\v,t)+G'i\v 1 t) 



2 
2 

w 



Re 
Re 



d d k 



(27r) d 



-i(k— k p )-r 



$* 2 



(2n) d 



3 i(k-k p )-r 



^ r (<ty>k(i)<¥k(0)) 



(37) 



where k = 2k p — k. In the presence of quantum fluctua- 
tions onl y 14 i 22 , the normal correlations are 



(<ty£(i)<tyk(0)) =27rft 2 73 2 |$| 



and the anomalous correlations an 



A u (cj) 



-3(5 



Wk(i)%(0)) = -27rft 2 7 .g$ 2 / du; 



„— tort C k 



(38) 



A k (w) 

(39) 

Here we have introduced the notation £k = £k + 2g|$| 2 — 
Hujp for the shifted single-particle dispersion of polaritons, 
k p and lo p respectively indicate the wave vector and the 
frequency of the pump. The denominator Ak is defined in 
terms of the eigenvalues ' of the Bogolubov matrix 

.9$ 2 



as 



, e k -ifvy/2 



ihj/2 



(40) 



AkH = {twj-E^){twj-E^ )(fiw-4 2) )(?^-4 2) ). 

(41) 

To better compare to the numerical results, we now re- 
strict our attention to one-dimensional case and we ex- 
tract analytic formulas for the most significant regimes 
considered in the previous subsection. 



1. Low- density parametric luminescence 

In the low-density limit |$| — > the Bogolubov eigen- 
modes tend to the single-particle energies with a negligi- 
ble blue-shift [Fig. [T](a, b)]. Comparing the expressions 
(|3H1) an( l (j3"9")l for respectively the normal and the anoma- 
lous correlations, it is immediate to see that the former is 
a factor g |$| 2 /fi,7 smaller and therefore negligible in this 
limit. By performing the frequency integral in (|37j) with 
the method of residuals, the spatial correlation function 
can be written as 



G {2) (r,t) ~ -47rgSgn(t)e- 7l * l/2 x 

e -ie k t/h 



Re I / dke l(k - k " )r 



e k + e k — 2hu> p — itry 



(42) 



In the parametric configuration (see Fig. Q] (a, b)), 
the dominating contribution to the integral over k is 
given by the wavevectors around the signal k s and the 
idler fcj = k s = 2k p — k s modes, where the normal and 
the ghost dispersions intersect and parametric processes 
are resonant. We can then approximate the integrand 
by replacing the energies by their first-order expansions 
around k = k s ; 



with group velocities 



(43) 



(44) 



In this way, we obtain 



G { 2 2) (r,t) ~ 87r 2 -f— Sgn{t)e-^ /2 Ssm(Kr - M): 



hAv 



x [6(tSR s )e- KR ° - 6(-tSRi(r, t))e Kfl ' (r ' t) ] , (45) 



(») 



where S = Sgn(Au s ) and Av g 

From this expression, it is immediate to see that the 
system of parallel fringes has a frequency Ml = e kg — hw v 
and a wavevector K — k s — k p determined by the in- 
terference of the signal/idler and the pump mode. The 
analytic form of the zero delay t = fringe pattern 
resulting from the combination of sin and 9 functions 
is responsible for the dip at r — that is visible in 
FigHIb). The temporal decay of the correlation occurs 
on the same time scale as the bare polariton decay rate 
7. The spatial decay away from the butterfly edges (we 



t) occurs on a length scale 



have set R s ^(r,t) = 1 
n- 1 = \Av g \h. 

Furthermore, the analytical formula (|4"5")) shows that 
no correlation is present inside the cone marked by the 
thin lines in the Fig. [2] and defined (for t > 0) by the 
condition R s < and Ri > (for t < the signs are 
exchanged). The physics behind this fact is illustrated 
by the simple geometric construction shown in Fig. |6] 
pairs of entangled signal/idler polaritons are generated 
at all times and positions by the parametric conversion 
of quantum fluctuations into real excitations. Signal and 
idler polaritons then propagate with group velocities re- 
spectively v^ 1 and transport the correlation to distant 
pairs of points. It is easy to see that the shaded region 
inside the cone can never be reached by such a process. 
An analogous reasoning was used in Ref<^ to explain the 
correlation signal observed in numerical calculations of 
Hawking radiation from acoustic black holes. As polari- 
tons decay at a rate 7, the same construction shows that 
the correlation signal has to decay in space with the char- 
acteristic length k^ 1 = Aiig/j. 

The perfect agreement between this analytic approxi- 
mation and the numerical result is highlighted in Fig. [2] 
(b. c). The lower polariton parameters g, 7, e k are cal- 
culated from the linear eigenmodes of the GP equations 
© and (fTU|) . As expected for a Bogolubov theory^ 4 -, for 
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r(nm) 

FIG. 6: Geometric construction to determine the pairs of 
points in space-time that are correlated in the low-density 
parametric luminescence regime. The shaded area highlights 
the region of space-time which is not correlated with the point 
(r = 0, t = 0). The point (r = 0, t = 0) is parametrically cor- 
related with the point a, as a result of the parametric emission 
event occurring at the point marked with a star. The same 
holds for the point /3. Solid lines mark the motion of signal 
(blue) and idler (red) particles. No parametric event can in- 
stead produce correlations between the point [r = 0, t = 0) 
and the point S located inside the shaded region. 



a given blue shift g |$| 2 the rescaled correlation function 
is proportional to the nonlinear coupling constant g. 

Before proceeding, it is interesting to note that for a 
pump at normal incidence k p = 0, one has k s = —ki and 
a vanishing relative frequency O = 0. The spatial fringes 
are then independent on time. From an experimental 
point of view, this configuration appears to be the most 
suitable one to observe the predicted features as it is less 
subject to the finite time resolution of photodetectors. 



2. OPO critical region 

For pump intensities just below the OPO threshold, 
fll 3" 1 2 ^ g\^OPo\ 2 = fry/2, the Bogolubov spectrum 
strongly differs from the single-particle one [Fig. [T] (e, 
f)]. In particular, the imaginary part of one of the two 
eigenvalues tends to zero in the vicinity*^ of fc = k s 



and k 



2k p — k s . In this case, it is easy to 



i one has 
(i) 1 



see that for wavevectors in the vicinity of fc, 

ME ( kl +k \ = ME ( "\ +k ] and Re[E^>_ k ] - 
In this region, the eigenvalues can be approximated with 
the expression 



(1,2) 



ih 

■±E +kVg(k-k s> i)- — [jTlo±li(k-k Sti ) 



E k k_».fc, i 

(46) 

where V g = d k E k \ k Jh = d k E k \ k Jh and fry = 2g |$| 2 < 
tvj, the equality holding exactly at threshold |$| 2 = 
\<&OPO \ 2 ■ Following the procedure already adopted in the 



previous subsection, we perform the frequency integra- 
tion in (|39p with the method of residuals and we expand 
the function in k around fc s ,;. Then, by retaining only 
the dominant contribution (which, as expected, diverges 
exactly at threshold 7 = 70) and using the relations (|46[) . 
we finally obtain the expressions 



2 i e -(7-7o)|*|/2 CQs{Kr _^ )x 

h h 



dk 



cos(kR) e 



-71 fc 2 I* 1/2 



7 - 7o + lik 2 



(47) 



and 



Airg 



x Re 




e -(7-7o)l*l/2 x 



e ifei?-7ife 2 |t|/2 

dk — x 

7 - 7o + 7i « 



AV g {s) k + ihj /2)e 
AV g (i) k + ih-f /2)e 



Kr+iE t/h , 
iKr-iEat/h \ 



with R = r- V g t, E s 4 = e ks , =F E , K 



(48) 
fc p , and 



0,0 



V a 



(s,i) 



This result can be further simplified by introducing the 



approximations E s ~ Ei ~ and AV & 



» 



AVc 



(i) 



0. 



In this case, the anomalous correlation disappears 



4irg 



-(7-70)1*1/2 cos (^ r _ Eot/h) x 



dk s'm(kR) 



-7ifc 1*1/2 



7 - 7o + lik 2 



= , (49) 



as the integrand is odd in fc, and the spatial correlations 
are dominated by the normal contribution ~ G 2 . 

— _____ 

The explicit expression of G\ is given in (|4"7|) : the 
Fourier transform of the product of a Gaussian and a 
Lorentzian function corresponds, in real space, to the 
convolution of an exponential and a gaussian function. 

At zero delay t = [Fig. [4] (b)], the Gaussian reduces 
to a delta-function in space and the correlation signal 
shows an exponential decay in space with a characteristic 
length 



7i 



7-7o 



(50) 



As expected on the basis of general arguments on phase 
transition, and previously observed in Monte Carlo simu- 
lations^, the characteristic length £ diverges as the OPO 
threshold is approached 70—5-7. In contrast to the ex- 
pression (|4"5j) for the weak intensity case, the zero delay 
t = fringe pattern now has a cos form with a simple 
bunching peak at r = 0. 

At longer times, the spatial width of the Gaussian 
grows as y/i, so that for short to intermediate distances, 
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the dependence is dominated by the Gaussian factor. 
This effect is clearly visible in Fig. [I] (c). The overall 
exponential decay in time occurs on a characteristic time 

t~-^—. (51) 

7-7o 

which again diverges as the critical point is approached. 

3. Non-parametric regime 

The non-parametric regime shown in Fig. Q] (g,h) 
corresponds to the case where the polariton and ghost 
branches do not intersect, i.e. + ^ for all k. 
This regime is generally realized when the frequency of 
the pump is very much red-detuned with respect to the 
renormalizcd polariton dispersion and, in our configura- 
tion, is fulfilled for pump intensities above the bistability 
loop. 

In this regime, the Bogolubov eigenmodes can be ap- 
proximated by the single particle dispersion blue-shifted 
by the interaction, i.e. E k — £ k — ifvy/2 and an equation 
formally equivalent to (|42p still holds. However, since 
the denominator remains finite for all k 7 no pole can be 
identified. However, the region of the polariton disper- 
sion that minimizes the denominator gives the dominant 
contribution. In the case of Fig. Q] (g), this happens at 
k = k p . It is therefore convenient to expand the energy 
dispersion at second order in k — k p . Performing this ap- 
proximation in the denominator, we obtain the following 
formula for the correlation function: 



G^(r,t) ~ ig(ir^- ~ 2)e~^/ 2 x (52) 

e ie k+kp t/H 



x Re < dke 



-ikr 



2e k + 2u p k 2 + iSgn(t)hj 



To further simplify this expression, we can expand also 
the energy exponents, which leads to 



G^(r,t) ~ 4 3 (^^- -2)e-tW 2 Re\e ii *p t / A x (53) 

e k p L 

piUpk 2 1 / h 



x / dke 



-ikR 



2e kp + 2u p k 2 + iSgn(t)/ry 



Jh and 



Here, we have set R = r — Vgt, v g p ^ = d k e k \ 
u p = 9 2 k e\ kp /2. 

The integral over k is of the Fresnel kind and describes 
the interference produced at the point (i?, t) by the differ- 
ent fc-modes with a gapped and quadratic dispersion. For 
zero time delay t = 0, the correlations have a typical anti- 
bunching character: they are everywhere negative and 
are strongest at r — 0. Further away, they monotonically 
tend to zero with an exponential law of characteristic 
length = y/ u p /e kp determined by the gap between 
the renormalizcd polariton and ghost branches 3 ^. 



The most apparent deviation between the analytical 
form (|5"3"|) and the numerical result shown Fig(5jc) con- 
sists of a tail in the analytic approximation that extends 
up to large distances. This has a simple interpretation: 
the quadratic approximation of the dispersion eliminates 
all bounds in the group velocity and predicts correlations 
at any distances. In contrast, the correct dispersion has 
an upper bound r>™ ax to the group velocity, which re- 
stricts the possible correlations to the |r| < v™ ax \ t\ re- 
gion marked by the thin lines in Fig. [5{a) and (c): this 
intepretation is confirmed by the much better agreement 
of the (|5"2"j) prediction where the group velocity is cor- 
rectly taken into account. 



IV. TWO-DIMENSIONAL AND DISORDERED 
SYSTEM 

In this final section we apply our model to the more 
general case of a two dimensional inhomogeneous system. 
In particular, we wish to investigate how the conclusions 
of the previous sections are affected by the presence of 
exciton and photon disorder. 

Realistic system parameters for a GaAs microcavity 
with N = 2 quantum wells are used, with Ml r = 3.5 
meV, h"f c = tvy x = tvy = 0.2 meV, zero exciton-photon 
detuning. For the nonlinear interaction constant we take 
g x = 1.5 x 10~ 3 meV /im 2 and g s = 0.5 x 10 -3 meV 
fim 2 . We consider a configuration where the pump is 
orthogonal to the cavity plane, k p = and we take Huj p — 
e c (0) = —3 meV. This orthogonal pump configuration 
is the most suitable one in view of experiments, as it is 
least affected by the temporal resolution of the photon 
detectors. 

The correlation pattern for different values of the pump 
intensity, corresponding to the low-density limit, the crit- 
ical OPO region and the non-parametric configuration 
are shown in Fig. [7] (a, b, c) for a two dimensional sys- 
tem in the absence of disorder. All the features discussed 
in the previous section for the ID case are still apparent. 
In particular, in the low-density case, the correlations 
disappear for distances smaller that |r mo a;(*)| = v g t, ±v g 
being the group velocities of polaritons in the signal and 
idler modes, and they decay exponentially. On the other 
hand, in the vicinity of the OPO threshold, correlations 
extend everywhere. In the non-parametric configuration, 
correlations are non vanishing only for distances smaller 
being the largest group velocity. 



J 9 



.max 

J g 



We assume white noise disorder for the exciton field 
of amplitude 2 meV and a gaussian-correlated disorder 
for the photon field, with amplitude V% is = 0.5 meV and 
a correlation length of £ c = 7 yum. The results of the 
numerical calculations are summarized in Fig. [7] (cl-i): 
the profile of the coherent photon field |$ c (r)| 2 in the 
three cases is shown in panels (g-i) and the correspond- 
ing spatio-temporal patterns of correlations are shown in 
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-20 20 -20 20 -20 20 

x (p.m) x (|im) x (|im) 



FIG. 7: Intensity correlation pattern G' 2 '(r,i) as a function 
of the spatial coordinate x and of the temporal delay t along a 
y — line for a 2D system. Panels (a-c) are for a uniform case, 
while panels (d-f) are for a disordered system. The pump in- 
tensity is varied from a low value (a, d), to a value close to 
the OPO threshold (b, e) and to a value well above the bista- 
bility loop (c, f) where parametric processes are forbidden. 
Panels (g-i) show the coherent photon field profiles |<l> c (r)| 2 
corresponding to the disordered case of panels (d-f). 

panel (d-f). The realization of the disorder potential is 
the same for the three values of the pump intensity. 

Even if it is responsible for large density modulations, 
the considered disorder is never able to destroy the near- 
field correlation pattern. Of course the pattern would 
eventually disappear if a much stronger disorder was con- 
sidered, that is able to fragment the coherent field in dis- 
connected parts. However, this latter situation is quite 
unusual for state-of-the-art GaAs microcavitics. 

From the comparison with the corresponding correla- 
tion patterns for a clean system [Fig. [7J (a-c)] , we can 
still appreciate a slight modification in the pattern. Dis- 
order is responsible for these modifications via three main 
effects: i) the exciton and photon resonances are broad- 
ened; ii) k-conservation is broken, which softens the con- 
dition for parametric processes; iii) the group velocity is 
no longer a well defined concept, so the contour of the 
butterfly shape in Fig. [7J (d) is smeared out with respect 
to panel (a). 

V. CONCLUSIONS 

We have developed a formalism to compute the second 
order spatial correlation function of polaritons in a pla- 
nar microcavity. This quantity directly transfer into the 
near-field intensity correlations of the emitted light. 

We have computed the spatio-temporal pattern of cor- 
relations in different pumping regimes, ranging from the 
parametric luminescence regime, to the critical region 
just below the parametric oscillation threshold, and to 
the strong pumping regime where a large blue-shift of 



the polariton modes is able to prevent parametric oscil- 
lation. 

For each regime, we have identified the key features 
of the correlation pattern. In the uniform case, we have 
compared our results to the predictions of approximate 
analytic models which provide a physical interpretation 
to the observed patterns. An orthogonal pump geometry 
appears as the optimal choice for the experimental ob- 
servation as it reduces the required temporal resolution 
of photo-detectors. We have verified that the correla- 
tion patterns are not qualitatively modified by a realistic 
disorder as long as the coherent polaritons remain dclo- 
calized in space. 

The conclusions of the present work confirm the ex- 
pectation that intensity correlations can be a very pow- 
erful tool to study the dynamics of quantum fields in 
condensed matter systems. A study of more complex 
geometries is presently under way. 
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VII. APPENDIX 

In this Appendix, we discuss the connection between 
the photon field inside the cavity and the emitted light. 
In particular, we show that the calculated in-cavity in- 
tensity correlations directly transfer to the near-field cor- 
relations of the emitted light. 

We consider a planar cavity along the x—y plane with a 
quantum well placed at z = and mirrors whose external 
surface is at z = ±z m . In Fourier space, the external 
field £" OM t(k,w), resulting from the transmission across 
the mirror is related to the in-cavity field Ei n (k,uj) at 
the quantum well position via the complex transmission 
coefficient r(k,w) = T l / 2 e 1 ^^: 

E out (r,z,t) = J " / ^( k .r-*VM*)(*-^) x 

x T(k,u})E in (k,oj) . 

For simplicity, we have neglected here the k- and in- 
dependence of the transmittivity T as we are resticting 
to frequencies well within the stop band of the mirror 
where T is almost constant (see FigJHJa)). 

For frequencies close to the k = cavity frequency ujq 
and at low wave vectors k, the phase of the transmis- 
sion coefficient can be accurately approximated by the 
expansion 

0(k,w) = 0(0, w ) + A t (w - u ) - ^-A z k 2 , (54) 
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with 



and 



x 10" 



z ~ c dk 2 



The transmitted field thus reduces to 
d 2 k 



(55) 



(56) 



Eout(r,z,t) ~ B 



(2n)< 



where we have used the relation (valid in the air) 
k z (k) = JujI/c 2 - k 2 ~ ujq/c - c/2uj a k 2 , 



(57) 
(58) 



and we have included the transmittivity T and a global 
phase into the multiplicative constant B. 

Without the mirror, the field E out (r, z,t) would be 
given by 



E out (r,z,t) 



d 2 k 



e lkr e" 



-zk' 



E m (k,t). (59) 



Comparison between Eqs. ([57]) and (|59j) shows that the 
only effect of the mirror on the external field is to give a 
time shift of A t and a longitudinal space shift A z . 

We now demonstrate that the expansion (|54l) holds 
for a typical high-quality GaAs microcavity with TV 
DBR's mirrors^. To this purpose, we employ the stan- 
dard transfer matrix approach^ to compute the k- and 
w-dependent complex transmission coefficient t(cj, k) 
across the top mirror placed between the cavity and the 
air. We consider the typical case where the DBR's mir- 
rors are composed by two alternate dielectric layers with 
index of refraction n\ and n 2 and length l\ and 1%, re- 
spectively, such as h/ni ~ £2/^2 — Ao/4, Ao = 2-kc/ujo 
being the cavity wavelength. The transfer matrix for the 
transmission from the cavity to the air is2£ 



M = TaaT^M&Tu 



(60) 



where T a2 , T\ 2 and T\ c are the matrices describing the 
transmission across the interfaces between the air and the 
dielectric 2, between the dielectric 1 and the dielectric 2 
and between the dielectric 1 and the cavity, while 



M12 = T 21 M 2 T 12 Mx , 



(61) 



is the transfer matrix describing the transmission across 
the periodic block formed by the two dielectric layers. 
Since 



M 



the transmission coefficient is given by 

_ det(M) 
Mm ' 



(62) 



(63) 




k (urn" 1 ) 



FIG. 8: Transmittivity (a) and transmission phase (b) across 
the A/4 DBR mirror between cavity and air as a function of 
the frequency u, for a given wavevector k = (solid line), 
k — 1 fim^ 1 (dashed), k = 1.5 /im -1 (dot-dashed) and k — 
2 /xm _1 (dotted). Parameters are h = 65 nm, h = 75 nm, 
ni = 3, n2=3.5, n c = 3.4, L c = 290 nm and N = 22. The 
cavity frequency ujo = nc/n c L c is located at the center of the 
mirror stop band, (c) Transmission phase across the same 
DBR mirror as a function of the wavevector k, for a give 
frequency oj — OJo = (thick solid line) and h(u — cjq) = ±2 
meV (thin lines). Dashed lines are obtained by fitting the 
numerical results with Eq. Q54p . which gives A 2 ~ 140 nm. 



Examples summarized in Fig. [5J In panel (a), we 
show the frequency dependent transmission amplitude 
T, for in-plane wave vectors ranging from k — and 
k = 2 /im" 1 . The relative variation is less than 10% in 
the interval of frequencies and wave vectors considered 
which validates the assumption underlying (). In panel 
(b) we display the phase of the transmission coefficient 
as a function of the frequency. For all values of the in- 
plane wavevector k, the linear dependence assumed in 
Eq. (f5"4"]) is accurately verified with the same A t ~ 5 fs. 
This result confirms the validity of the expansion f|54[) . 
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In panel (c) , we show the dependence of the phase <f> on 
the wavevector k, for different values of the frequency 
hui = —2, 0, 2 mcV. Again the dependence assumed in 
Eq. (|54p is verified with the same value of A 2 ~ 140 nm. 

It is important to note that this small spatial shift 
A 2 does not scale with the number N of layers in the 
DBR mirror, as it could be intuitively supposed. This 
result originates from the fact that within the stop band 
the transfer matrix Mi 2 (and consequently M/^) has real 
eigenvalues, and thus does not induce any extensive phase 



shift. 

The results of this appendix show that the only effect 
of the mirror on the near-field correlation pattern is the 
following: the light detected outside of the cavity ap- 
pears to be generated at a slightly earlier time and at 
a slightly displaced longitudinal position as compared to 
the quantum well position. In the experiments one has 
therefore simply to focus the optical detection on this 
shifted position. 
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